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Abstract 

Central moment lattice Boltzmann method (LBM) is one of the more recent developments 
among the lattice kinetic schemes for computational fluid dynamics. A key element in this 
approach is the use of central moments to specify collision process and forcing, and thereby 
naturally maintaining Galilean invariance, an important characteristic of fluid flows. When the 
different central moments are relaxed at different rates like in a standard multiple relaxation time 
(MRT) formulation based on raw moments, it is endowed with a number of desirable physical 
and numerical features. Since the collision operator exhibits a cascaded structure, this approach 
is also known as the cascaded LBM. While the cascaded LBM has been developed sometime ago, 
a systematic study of its numerical properties, such as accuracy, grid convergence and stability 
for well defined canonical problems is lacking and the present work is intended to fulfill this need. 
We perform a quantitative study of the performance of the cascaded LBM for a set of benchmark 
problems of differing complexity, viz., Poiseuille flow, decaying Taylor- Green vortex flow and 
lid-driven cavity flow. We first establish its grid convergence and demonstrate second order 
accuracy under diffusive scaling for both the velocity field and its derivatives, i.e. components of 
the strain rate tensor, as well. The method is shown to quantitatively reproduce steady/unsteady 
analytical solutions or other numerical results with excellent accuracy. The cascaded MRT LBM 
based on central moments is found to be of similar accuracy when compared with the standard 
MRT LBM based on raw moments, when detailed comparison of the flow fields are made, 
with both well reproducing even small scale vortical features. Numerical experiments further 
demonstrate that the central moment MRT LBM results in significant stability improvements 
when compared with certain existing collision models at moderate additional computational cost. 
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I. INTRODUCTION 



Early developments in the area of computational fluid dynamics (CFD) have focused on 
the solution of the classical discretizations of the continuum description of fluid motion. 
During the last two decades, there has been much interest and effort in the development of 
schemes that derive their basis on a more smaller scale picture involving particle motion, 
which may be classified as mesoscopic methods. One of the most promising of such ap- 
proaches is the lattice Boltzmann method (LBM) [1-3J. Based on kinetic theory, it involves 
the solution of the lattice Boltzmann equation (LBE), which specifies the evolution of the 
particle populations along discrete directions, which comprise the lattice. This evolution 
involves a Lagrangian free streaming process along such lattice links and a local collision 
step specified as a relaxation process. Various elements involved in these two simple steps 
are constructed based on symmetry considerations, while obeying certain conservation con- 
straints, in such a way that they recover the dynamics of fluid flow in the near incompressible 
limit. The resulting scheme has a number of desirable features. These include the ability 
to naturally represent complex fluid physics such as multiphase and multicomponent flows 
based on kinetic theory, amenability to parallelization due to the locality of the method and 
representation of flow through complex geometries. Furthermore, due to the exact conser- 
vation in the streaming step and machine round-off conservation in the collision process, it 
has considerably low numerical dissipation for a second-order numerical scheme Due to 
such competitive advantages, the LBM has found applications in the simulation of a wide 
range of fluid flow problems [lH3]. 

Since the LBM is usually developed by means of a bottom-up strategy, there is certain 
level of flexibility in the construction of its various elements to recover the macroscopic fluid 
motion. In particular, the choice of a suitable collision model can have profound influence 
on the fidelity as well as the stability of the approach. As such, the construction of the 
collision step has been the subject of considerable attention since the inception of the LBM. 
The simplest among these is the so-called single-relaxation-time (SRT) model [HIE], which is 
based on the Bhatnagar-Gross-Krook (BGK) approximation [?]. While it is popular, it has 
limitations in the representation of certain flow problems and is generally prone to numerical 
instability, particularly at high Reynolds numbers. A major development to address these 
aspects is the moment approach pj , which has been constructed based on multiple relaxation 
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times (MRT) in particular to significantly improve the numerical stability [9]. While it is 
related to its precursor involving a more general relaxation approximation piU [TT] . the 
characteristic difference being that it performs collision in an orthogonal moment space 
leading to an efficient and flexible numerical scheme. This moment approach, which is 
designated as the standard MRT formulation in this paper, has recently been studied and 
compared with some of the other collision models in detail [1^. A simpler version that is 
intermediate between the SRT and MRT model is the so-called two-relaxation-time (TRT) 
model [13], in which the moments of even and odd orders are relaxed to their equilibrium 
at different rates. This, along with the MRT model, can be adjusted such that it results 
in a minimization of undesirable discrete kinetic effects near walls. Another significant 
development is the so-called entropic LBM [Hj. It involves an equilibria, which is based 
on a constrained minimization of a Lyapunov-type functional. By modulating the collision 
process through enforcing entropy involution locally, this approach aims to maintain non- 
linear stability. This approach has resulted in a number of simplified variants recently 
[IS]. 

An important physical feature of the fluid motion is that their description be independent 
of any inertial frame of reference (e.g. [IZ])- This invariance property, which is termed 
as the Galilean invariance, should be satisfied by any model or numerical scheme for its 
general applicability. Furthermore, it has recently been shown that stabilization of classical 
schemes for compressible flow can be achieved when they are specifically constructed to 
respect this physical property [TSH2U] . Keeping these general notions in mind, Galilean 
invariance can be naturally prescribed in the LBM when its various elements are represented 
in terms of the central moments, i.e. moments obtained by shifting the particle velocity by 
the local fluid velocity. That is, any dynamical changes due to the collision process and 
impressed forces can be represented in terms of suitable variations of a set of such central 
moments. In particular, a collision model based on the relaxation of central moments was 
constructed recently [21]. The model exhibits a cascaded structure, which was later shown 
to be equivalent to considering a generalized equilibrium in the lattice or rest frame of 
reference [22]. These central moments can be relaxed at different rates during collision 
leading to a cascaded MRT or central moment MRT formulation, whereas by contrast the 
standard MRT formulation considers raw moments. A systematic derivation of this approach 
by including the effect of impressed forces based on central moments was presented in [23] . 
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This leads to considering generalized sources, analogous to the generalized equilibrium in 
the rest frame of reference. They also presented a detailed Chapman- Enskog analysis of the 
cascaded MRT LBM for its consistency with the macroscopic fluid dynamical equations of 
motion. This approach was further extended to various lattice models in three-dimensions 
in [2lj, in the cylindrical coordinate system for axisymmetric flows in [25] and for accounting 
of non-equilibrium effects in [26] . 

Prior work on the cascaded LBM as discussed above have focused mainly on method 
developments or their mathematical analysis, with little attention towards their numerics 
except for few validation cases. In particular, a detailed numerical study of the properties of 
the cascaded LBM for established benchmark problems and also their performance against 
other LBM approaches is lacking. The focus of the present work is intended to fill this 
gap by presenting a systematic study of the numerical properties of the cascaded LBM, 
viz., grid convergence, accuracy and stability for various canonical problems of differing 
complexity in terms of flow features and temporal evolution. Establishing the reliability 
and merits of the method in quantitative terms could provide confidence in their extension 
and applications to various complex flow problems of interest. To study the numerics of 
the cascaded LBM, we consider the Poiseuille flow, decaying Taylor-Green vortex flow, and 
lid-driven cavity flow, for which either analytical solutions or detailed prior numerical results 
are available for comparison. Much of the literature on the LBM with other collision models 
on grid convergence studies have focused only on those for the velocity fleld. In this work, 
we present numerical results on the grid convergence of the cascaded LBM for the velocity 
fleld as well as its derivatives, i.e. the strain rate tensor. Furthermore, an advantage of the 
kinetic schemes such as the LBM is that the strain rate tensor can be computed locally in 
terms of non-equilibrium moments. In this work, we also present a direct comparison of 
the results obtained using the non-equilibrium moments of the cascaded LBM with those 
involving the flnite differencing of the velocity fleld at various locations for the lid-driven 
cavity flow problem to assess their quantitative accuracy. It may be noted that a detailed 
comparison study of the SRT and the standard MRT models have recently been performed 
in Thus, in this work, we present a quantitative accuracy comparison between the 

standard MRT LBM and the cascaded or central moment MRT LBM for the lid-driven 
cavity flow. Finally, we will discuss the numerical stability performance of the various LBM 
schemes for the above benchmark problem. 
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The paper is organized as follows. Section |TT] presents the details of the particular version 
of the cascaded MRT LBM used in this work. In Sec. |III[ the results of the grid convergence 
study of the cascaded MRT LBM together with the raw moment based standard MRT LBM 
for the three benchmark problems are discussed. Subsequently, the quantitative accuracy of 
the cascaded LBM is demonstrated by making detailed comparison with either analytical or 



other numerical solutions for the above problems in Sec. IV In Sec. |V| numerical stability 
test results are presented for the lid-driven cavity flow using the SRT LBM, standard MRT 



LBM and cascaded MRT LBM. Summary and conclusions of this work are given in Sec. VI 



II. CASCADED LATTICE BOLTZMANN METHOD 

We will now discuss the main features of the cascaded LBM. Similar to the standard MRT 
LBM, the cascaded MRT LBM also performs collisions in moment space, but these moments 
are obtained by shifting the particle velocity by the local fluid velocity, i.e. using central 
moments. As a result, the approach can naturally maintain Galilean invariance. Central 
moment relaxation process was specified in [21], which was re-interpreted by considering 
generalized equilibrium in [22]. Its detailed mathematical consistency analysis in a MRT 
formulation with forcing was carried out in [23] • The computations of the cascaded LBM 
are actually performed after transforming the central moments into raw moments by means 
of a binomial formula. In this work, the specific formulation of the cascaded LBM given 
in [23], whose details are somewhat different from that given in [21], is used. This is briefly 
discussed in what follows. 

In this work, the standard two-dimensional, nine velocity (D2Q9) lattice is employed. 
We consider the usual bra-ket notations in the description of the method as it provides a 
convenient representation. That is, we consider the depiction of vectors as (0| and where 
(01 represents a row vector of of any state in the corresponding direction (0o, 0i, 02, ■ ' ' , ^s) 
and 10) represents a column vector (0o, 0i, 02, ■ ■ ■ ,4>%Y- The inner product Y^^a=Q't>a'^a is 
then denoted by {(t)\^)- As the cascaded LBM is a moment approach, we need a set of nine 
linearly independent moment basis vectors for its specification. The (raw) moments of the 
distribution function of different orders can be defined as Y^a=o^^x^^yfcf Here, a is the 
discrete particle direction, and m and n are integers. Thus, a set of nine linearly independent 
nonorthogonal basis vectors obtained using the monomials e™^.e" in an ascending order can 
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be written as 



|p) = ||e-;|°) = (1,1, 1,1, 1,1, 1,1,1)^, 
|e„,) = (0,1, 0,-1, 0,1, -1,-1, If, 
|e„,) = (0,0,1,0, -1,1,1,-1, -If, 
|eL + e^,) = (0,l,l,l,l,2,2,2,2f, 

\eL-ely) = (0,1, -1,1, -1,0, 0,0, Of 
Kxtay) = (0, 0, 0, 0,0, 1,-1,1, -if , 
|eLeay) = (0,0,0,0,0,l,l,-l,-lf, 
|eaxe^J = (0, 0,0, 0,0, 1,-1,-1, If, 



(1) 



|eLO = (0,0,0,0,0,l,l,l,l)^ 
This can be transformed by means of the Gram-Schmidt procedure into an equivalent set 
of orthogonal basis vectors, which provides a computationally more efficient and convenient 
setting for the description of the method. As a result, we have the following orthogonal 
set ESI: 



K2) 



K4) 



\P), 

I Coa;) ) 



^3)=3|eL + e^,)-4|p), 



I 2 2 \ 

I ^ax^ay") 1 

'^\^ax^OLy) "I" 2|eQ,j^), 

'^\^ax(iay') ~^ '^\^ax^ 1 



(2) 



^8)=9|eLe^,)-6|eL + e^,)+4|p). 

Collecting the above set of vectors as a matrix /C, it immediately follows that /C/C^ is 
a diagonal matrix, owing to orthogonality. This orthogonal matrix /C can be written in 
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(3) 



component form as 

/C = [\Ko), \K,), \K,), \Ks), \K,), \K,), \K,), \Kj), \Ks))] 

100 -4 00004 
11 0-11 2-2 

10 1-1-10 2 0-2 
1-10-11 0-2-2 

10- 1-1-10-20 -2 

11 1 2 1-1-11 
1-11 2 0-1-11 1 
1-1-12 1 1 1 1 

11- 12 0-11-11. 

To specify the collision step and forcing, we need the central moments of the local 
equilibrium and sources, which can be obtained as follows. First, the local Maxwell- 
Boltzmann distribution function in continuous particle velocity space {^x,^y) is written as 



2^exp 



2c2 



where is the speed of sound. Typically, 



1/3. Based on this, the continuous central moments of the equilibrium of order (m + n) 



can be defined as H^nyn 



IZo IZo /^(^- - ^^)"'i^y - ^vTdi^diy, which yields 



in 



M 



\ _ (Y\M Y\M Y\M Y\M Y\M TT, 

(p,0,0,c^p,c^p,0,0,0,ctpf. 



xyi xxyi xyyi xxyy 



(4) 



Considering that the impressed forces only influence the fluid momentum, the central mo- 
ments of the sources of order (m + n) due to a force field (F^., Fy) defined by T-^myn = 
I^oc I^oo — Ux)"^{^y — UyYd^xd^y^ where A/-^ is the change in the distribution func- 



tion due to force fields, can be simply written as [23 

pj- \ ^ pj- pj- pj- fr \T 

^myU I (^J- Q , i ^ , J- y , i J. yy^ J. ^y^ J. J, J. , J- ^ 1/ J/ ' XXyy/ ' 



(5) 



= (0,F„F„0,0,0,0,0,0)^. 

Based on the above continuous central moments, the elements of the cascaded LBE can 
be formulated. Using the trepezoidal rule representation of the source term, the cascaded 
LBE can be written as [231 



+ eJt^t + 5t) = fa{x, t) + ^lii^g^t) + ^ [Sa(x,t) + 5'c,(£+e-,,t+5,)] . 



(6) 
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Here, the collision term Vf^ can be represented as Qp^ = 0^(f, g) = (/C-g)^, where f = \ fa) = 
(/o)/i)""" ) /s)^ is the vector of distribution functions and g = \ga) = {go,di^''' ^ds)^ 
is the vector of unknown collision kernel to be obtained later. Owing to the cascaded 
nature of the central moment based approach, it satisfies the following functional relation 
'ga = ^a(f , ^^), /3 = 0, 1, ■ ■ ■ , a — 1. The discrete form of the source term Sa in the cascaded 
LBE given above represents the influence of the force field {F^, Fy) in the velocity space and 
is defined as S = \Sa) = {Sq, Si, S2, ■ ■ ■ ,Ss)'^. Noting that Eq. ^ is semi- implicit, by 
using the standard variable transformation / = /« — ^•S'^, its implicitness can be effectively 
removed. This yields 

7„(x + eJt,t + 6t) = Ja{x,t) + fi^(^^i) + Sa(x,t)- (7) 

The derivation of the collision term, i.e. the collision kernel g and the source term S 
involves matching the discrete central moments and the continuous central moments of 
equilibria and sources, which are specified above, of all orders supported by the lattice set. 
We designate this step as the Galilean invariance matching principle. First, the discrete 
central moments of the distribution functions and sources of order (m + n) can be defined, 
respectively, asK^^^yn = {{eax-Ux)"^{eay-Uy)''\fa) and a^myn = {{eax-Ux)"'{eay-Uy)"-\Sa) ■ 
Also, in terms of the transformed distribution functions we define K^myn = {{eax — Ux)^{eay — 
Uy)"'\f^), which satisfies and similarly for the local equilibria 

Ti^myn = {{cax — Ux)^{eay — Uy)"\fa )• Then, the Galilean invariance matching principle 
reads 

n,^myn — ll^niyTi, J 
Cx^y" = T^myn. (9) 

This immediately specifies the various discrete central moments. Hence, we get 

l^reg \ _ /-^eq -^eq -^eq -^eq -^eq -^eq -^eq -^eq -^eq \T 



xxyy J 



and 



= (p,0,0,c^p,c^p,0,0,0,c^pf , 

\(7^myn^ ('^0) ^x^ ^yi ^xxi ^yyi ^xyi ^xxyi ^xyyi ^xxyy) 

= (0,F,,F„ 0,0, 0,0,0, Of, 



\ 1—^1 —^1 —^1 —^1 —'^Q ^eq ^e<7 

iKj-myn/ [Kq , , K,y , K,xx^ ^yy ^xy^ ^xxy^ ^xyyi ^xxyy) ' 

= {p,-IFx, -^Fy, cIp, cIp, 0, 0, 0, ctpf. 



(10) 



(11) 



(12) 



The next important step is to transform all the above discrete central moments in 
terms of raw moments, which can be readily accomplished by means of the following bi- 
nomial formula: ((e„. - «J-(e„^ - «^)"|(^) = (e™ + (e™ [E^li C;e^7(-l)X] Iv') + 

(e^[E™iCreari-l)X]lv^> + ([E::iQ"<r(-l)X] [Y.UC]e^:^yK-^yui]W), where 
Cq = pi/ ~ ?)0- Thus, we obtain the following discrete raw moments of sources a'^myn 
as 



^0 = 


{Sa\p) = 0, 


K = 






('S'al^-aj/) — 




- i^al^ax) ~ 


Ky = 


- {Sa\G^y) — 


^xy ~ 




^xxy 


— i^al^ax^ocy 


^xyy 


{^al^ax^ay 



yi 

IF'x'^X-i 



(13) 



^y^x ~^ '^■P'x^x^yi 
■^x^y ~l~ '^Fy^y^xi 

'-^xxyy ~ {^oi\^ax^ay) ~ '^■^x'^x'^y + ^-FyMyH^,. 

Based on the above, we now obtain the source terms projected to the orthogonal moment 
basis vectors, i.e. {K^\Sa), /? = 0, 1, 2, . . . , 8. This intermediate step is needed to obtain the 
source terms in the velocity space. It immediately follows that 







= 0, 


m\ = 




— F 

— ^ X-i 


ml = 


{K2\S^) 


J- yi 


ml, = 




= &{FxUx^ FyUy), 


ml = 


{Ki\S^) 


— '^{FxU'x ■^y'^y)i 


ml = 


{K,\S^) 






{K^\Sa) 


= (2 - 'iul.)Fy - dFxUxUy, 


my = 


{K7\Sa) 


= (2 - 3ul)Fx - QFyUyUx, 


^8 = 


{K^\Sa) 


= 6 [{3ul - 2)FxUx + {3ul - 2)FyUy] 



Equivalently, this can be written in matrix form as )C^S 
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(/C • S)„ 



{{Ko\Sa) , {Ki\Sa) , {K2\Sa) , . . . , {Ks\Sa)) = (mg, mf , m|, . . . , m|)^ = |m^). By exploit- 
ing the orthogonal property of /C, i.e. /C~^ = K,^ ■ D~^, where the diagonal matrix is 
D = diag((ii'o|-^o) , {Ki\Ki) , {K2\K2) , . . . , {Ks\Kg)), we exactly invert the above to obtain 
the source terms in velocity space Sa as 

1 



So 
Si 
S2 
S3 
S4 
S5 
Se 
Sr 
S. 



9 
1 

36 
1 

36 
1 

^36 
1 

36 
1 

36 
1 

36 
1 

36 
1 

'36 



{-ml + ml), 



{6m{ - + 9ml + ^1^7 - Smg), 

(6m2 - ^3 - 9ml + 6"^6 ~ ^mg) , 

(-6m{ -ml + 9ml - - 2m^^) , 

(-6m2 - ml - 9ml ~ ^'^e ~ ^mg) , 

{6ml + ^1^2 + 2?7i3 + - Sm^ - 3m^ + m^) , 

[-6m{ + 6m2 + 2ml ~ ~ ^'^e + '^''^7 + "^s) ' 

(-6mi - + 2ml + 9^5 + + Sm^ + m^) , 

{6ml - 6^2 + - 9ml + - 3m7 + ml) . 



(14) 



The discrete raw moments of the transformed distribution functions K^myn, which will be 
needed in the evaluation of the collision kernel, can be conveniently written as follows: 





= (/JP) 


= P, 




ifal^ax) 


— P^x 2^^' 




ifal^ay) 


= puy — 




if le^ ) 

\J a\ axl 


/ {1.3,5.6,7,8} 

-[ ? 




if ) 

\J a\ ayl 


/ {2,4,5,6,7,8} 



(15) 



f at 
fai 
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'{5,7} {6,8} 
^xy = (faKxeay) = \ J2 ~ ^ \ ^ ^ 

'{5,6} {7,8}^ 
^xxy ~ if al^ax^ay) = j ^ ] ~ ^ ] j ® fay 



{5,8} {6,7}^ 




^xxyy \J al^ax^ay/ 

where we have used (a Ea +^ + " " " ) ®7a = «(7ai +7^2 +7^3 + " " " ) +^(7/3i +7/32 +7/33 + 
•■■) + ■■■, with A = {«!, as, ■ ■ ■}, B = {(3i, (32, f^s, ■■■},■■■ , as a compact summation 
operator for ease of presentation. Furthermore, the raw moments of the colhsion kernels 
Eq,(^ ■ s)aGax^ay = J2/3{^i3\^ax^ay)di3 cire needed in its construction. Colhsion invariants 
of conserved moments imply = = ^2 = 0. Exploiting the orthogonal property of the 
matrix JC, the non-conserved moments of at higher orders, i.e. (3 = 3, 4, ■ ■ ■ ,8 can be 
obtained as follows [23]: 

J](/C-g). = ^(ir/3|p)?/3 = 0, 

a p 

^(/C ■ %)aeax = ^ (Kpleax) 9^ = 0, 
a p 

^(/C ■ g)aeay = {Kf}\eay) 9/3 = 0, 

a /3 

■ g)aeL = (^f^\^lx) 9P = 6?3 + 2^4, 
^(/C ■ %Uly = J2 (Kpl^ly) h = 6?3 - 2?4, (16) 



J5, 



a /3 

^(/C ■ gjaCaxeay = ^ {Kf}\eaxeay) 9/3 = 4^5 
a 13 

XI ■ g)"eLeay = (^/^I^Leaj/) ?/3 = -4?6, 

a 13 

XI • g)«eaxe^y = X (^/slea^e^y) ?/3 = -4^7, 

a p 

Using the above, the collision kernel gp of the cascaded collision operator = fi^(f , g) 
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{}C-g)a can be obtained as follows. Starting from the lowest order central moments that are 
non-coUisional invariants (i.e. K^x and higher), they are successively set equal to their local 
attractors based on the transformed equilibria. This step provides tentative expressions for 
based on the equilibrium assumption. This is then modified to allow for relaxation process 
during collision. That is, they are multiplied with corresponding relaxation parameters |2T] . 
In this step, care needs to be exercised to multiply the relaxation parameters only with those 
terms that are not yet in post-collision states (i.e. terms not involving /3 = 0, 1, 2, . . . , a — 
1) for a given See [23] for various details involved in this procedure. Here, we summarize 
the final expressions of the non-conserved collision kernels, which are given as follows: 



?4 = ^ I P{UI - - (k^^ - 2yy) - ^(ct^.^ - d'yy) ^ , (iS) 



?5 = ^ I pUxUy - K^y - ]^d',^y \ , (19) 



^ We f 2 -' -' -' 1 - 1 1 / - \ 

9& = ^\ '^P^xUy + f^xxy - '^U^'^xy " UyK^^ - -(T^xj/ > " 2^2/(3^3 + 9i) 

-2«.?5, (20) 

^ 1^7 f 2 - - - 1 1 / \ 

9t = ^{ '^PUxUy + K^yy - 2UyK^y " ^ ^ ^ " "^^2/^ ^ - -U^^^Q^ - Qi) 

-2uyg5, (21) 
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'^xxyy '^'^xl^xyy '^'^y'^xxy ~^ ^x'^yy ~^ '^y'^xx 



xllyl^ 



y'^xy 



)fxxyy } - 2?3 - ^My(3?3 + 9^ " ^Mx(3?3 " 9^) 



-4uxUyg5 - 2uyge - 2m^^7. (22) 

In the above, where /3 = 3, 4, 5, . . . , 8, are the relaxation parameters, satisfying the usual 
bounds < < 2. When a Chapman-Enskog expansion [27] is applied to the cascaded 
LBM, it can be shown to recover the Navier-Stokes equations with the relaxation parameters 
Us = and W4 = cjs = u'^ controlling the fbulk and shear viscosities, respectively (e.g., 
^ ~ '^s (w^^ ~ ^)) 123] ■ The rest of the parameters can be adjusted independently improve 
numerical stability. In this work, = 1^5 = ^ is selected based on the specified kinematic 
viscosity, while the rest of the relaxation parameters are set to 1. 

The cascaded LBE can now be re-written in the form of the usual stream-and-collide 
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procedure, leading to the following two steps: 

7a t) = 7„ (f , t) + + S^^,,t) , (23) 

J^{x + ea,t + 6t)=J^{x,t), (24) 

where the symbol "tilde" (~) in the above equations refers to the post-collision state of the 
distribution function. Expanding the collision term in the first step, the components of the 
post-collision distribution function can be explicitly written as 

70 = 7o + [do - 4(?3 - 9s)] + So, 

71 = 7i + [do + 91-93 + 91 + 2(?7 - gs)] + Si, 

12 =l2 + \9o + 92-93-94 + 2(?6 - ds)] + S2, 

13 = 73 + [90 -91-93 + 94- 2(^7 + gs)] + S3, 

1a = 1a + \go - ?2 - ?3 - ?4 - 2(^6 + gs)] + Sa, (25) 

75 = 75 + [?0 + ?1 + ?2 + 2^3 + ?5 - ?6 - ?7 + ^s] + S^, 

Jg = 7e + [go -gi + g2 + 2^3 - ?5 - ?6 + ?7 + gs] + Sq, 

Jy =77+ [go - ?1 - ?2 + 2^3 + ?5 + ?6 + ?7 + ds] + S7, 
Js = 7s + [go + ?1 - ?2 + 2^3 - ?5 + ?6 - ?7 + gs] + Ss- 

The hydrodynamic fields, i.e. the fluid density and the velocity then follow from taking the 
zeroth and first moments of the distribution function, yielding 

8 

p = E7a = (7Jp), (26) 

a=0 

8 _ I _ I 

pUi = J^/aCai + -Fi = ifjeai) + -Fi,i = x,y, (27) 

a=0 

and the pressure p satisfies p = c^p. A particularly useful feature of kinetic schemes such as 
the cascaded LBM is that the strain-rate tensor can be computed locally from a knowledge of 
the non-equilibrium moments. In fact, this can be shown by means of the Chapman-Enskog 
analysis, which was performed on the cascaded LBE in [23]. Setting the components of the 
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momentum as jx = pUx and jy = puy, such an analysis shows [23] 



ft 



(neq) 



(neq) 



(neq) 



3^3 
2 

1 



{dxjx + dyiy) , 

{dxjy + dyjx) , 



(28) 
(29) 
(30) 



where fj^'^'^^ 



— fp'^ are the non-equihbrium raw moments. Specifically, /a 



^yy, and % = K^y, whose equilibria are f^'^ = 2/3p + p{ul + ul), Jl^ = p{ul-ul), 



and f^'' = pUxUy, respectively [23]. It thus follows that 



9xjx 
9xjy ~l~ dyjx 



3^3 

2 

3a;4 

3^5 



8 



'-a=0 



P + PUx 



'-Q = 



o^ax^ay pUxUy 



(31) 
(32) 
(33) 



These specific expressions will be exploited in the numerical study of the cascaded LBM in 
the remainder of this paper. In the sections that follow, we will present the results obtained 
with the cascaded LBM for a set of benchmark problems to assess its numerical properties 
in terms of grid convergence, accuracy and stability. 



III. GRID CONVERGENCE STUDY ON THE BENCHMARK PROBLEMS 

We first perform a numerical study involving grid convergence for canonical flows includ- 
ing a steady 2D Poiseuille flow, a time-dependent 2D decaying Taylor-Green vortex flow, 
and a 2D lid-driven cavity flow characterized by various complex features. In the various 
figures presented in this section, the symbols represent the computed solution using the 
cascaded MRT LBM, the thin solid lines are the resulting slopes representing changes in the 
relative errors as the grid resolution increases, and the thick solid lines are the ideal slopes 
corresponding to second-order accuracy. In this work, a diffusive scaling is applied to per- 
form the convergence tests [28]. According to this scaling, the errors due to compressibility 
effects decrease at the same rate as the errors due to grid discretization thus prescribing a 
consistent limit process to represent incompressible flow. That is, the velocity scales in the 
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same proportion as the length scales. Equivalently, this means that the ratio of the Mach 
number and the grid Knudsen number remains constant for different grid resolutions, i.e. 
Ma/Kn = constant. 

A. 2D Poiseuille Flow 

The 2D Poiseuille flow is first considered. The flow is between two parallel plates of infinite 
length in the streamwise direction subjected to a constant body force. A periodic boundary 
condition is applied at the inlet and the outlet and a no-slip boundary condition at the solid 
boundaries by employing the standard half-way bounce back approach. The grid convergence 
is established by considering the following resolutions consisting of 3 x 24, 3 x 36, . . . , 3 x 192 
lattice nodes under diffusive scaling. The relaxation time for shear modes is set to r = 0.55 
that specifies and u^. The rest of relaxation parameters are set to unity. The flow is 
driven by a constant body force with the components specified to yield desired condition 
(see below) and Fy = 0. This classical flow problem has the well known parabolic profile 
as the analytical solution given by u{y) = Umaxi^ — y'^/L^), where Umax = is the 

maximum velocity occurring midway between the plates, u is the kinematic viscosity related 
the to relaxation time r as given in the previous section, and L denotes the half-width 
between the plates. Figure [T] illustrates the relative global errors between the computed 
solutions obtained using the cascaded MRT LBM and the analytical solutions for such flow 
at different Reynolds numbers of 100, 200 and 400. The relative global error, which quantifies 
the difference between the computed and analytical solutions, is defined as 

Relative Error = ' ""-M ^ (34) 

where Uc,i and Ua,i are the computed and the analytical solutions, respectively, and a standard 
Euclidean norm is used in the above measurements. It is seen that the relative errors have 
slopes of almost equal to 2.00, which tells that the cascaded MRT LBM is well-posed second- 
order accurate for this problem. In addition, the relative errors are seen to slightly increases 
as the Reynolds number increases. 
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Grid Resolution 



FIG. 1: Grid convergence of the cascaded MRT LBM for the vefocity field in a 2D Poiseuille flow 
with constant body force under diff'usive scahng. 



B. 2D Decaying Taylor-Green Vortex Flow 

The second problem considered is the decaying Taylor-Green vortex [22], which is a 
2D unsteady flow induced by a prescribed initial vortex distribution and decaying due to 
fluid viscosity. The fluid domain is a square of side 27r with no inflow/outflow and wall 
boundaries. The initial condition is set to be periodic array of vortices in both x and y 
directions as follows 



u{x, y,0) = — uq cos{kx) sm{ky), 
v{x, y,0) = + Uq sin(A;x) cos{ky), 



P{x, y, 0) =Po 



u 



. — cos {2kx) + cos{2ky)) 



(35) 
(36) 

(37) 



where k 



271 

N 



is the wavenumber, uq and po are the initial values for velocity and pressure, 



respectively. Here, N is the number of grid nodes in each direction. The temporal evolution 
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has the characteristic time scale given by T = Since there is no external energy supplied 
and because of the presence of fluid viscosity, the velocity field will decay with time due to 
fluid viscous dissipation. There exists an analytical solution for this problem which is a 
solution of the Navier-Stokes equations in a periodic domain and given by 

u{x, y,t) = — uq cos{kx) sin(%)e~^'^^'^*, (38) 
v{x, y,t) = + uo sin(A;x) cos(A;?/)e~^^ (39) 



p{x, y, t) =po - ^ 



cos(2/cx) + cos(2%) 



e-'""''. (40) 



Furthermore, the components of the strain rate tensor also satisfy the following explicit 
analytical solution: 

Sxx=^ = kuo sm{kx) sm{ky)e-^''^^^ (41) 



du 
dy 

1 f du dv 



Syy — Q_ — Sxx (^^2) 



In this test, the Reynolds number of the flow is set to Re = ^ = 14.4, where / = 2% 
is the length of the domain. A periodic boundary condition is applied to all the sides of 
the domain. We consider the following parameters in our grid convergence study: r = 0.55, 
k = 1,2 and Uq = 0.01. Applying the diffusive scaling, we obtain the relative global errors 
between the computed and the analytical solutions for the grid resolutions of 24 x 24, 48 x 48, 
96 X 96, 192 X 192 for a representative time t = 30. IT. In Fig. [2]shown are the relative errors 
for the u-velocity component, which have the slopes of 1.99 and 1.98 for the wavenumbers 
k = 1 and k = 2, respectively. Figure [3] shows the relative errors for the only independent 
strain rate tensor component Sxx with the slopes of 1.99 and 1.98 as well for the above two 
wavenumbers. Thus, it is evident that the cascaded MRT LBM is second-order accurate 
not only for the velocity field, but also for the components of the strain rates as well. This 
finding is consistent with a recent study with the SRT LBM for this problem [30] . 

C. 2D Lid-driven Cavity Flow 

Finally, the 2D lid-driven cavity flow is considered, whose geometric simplicity is con- 
trasted by various complex flow features. It is generally considered a standard benchmark 
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FIG. 2: Grid convergence of the cascaded MRT LBM for the velocity field in a 2D Taylor-Green 
vortex flow with k = 1 and k = 2. 
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FIG. 3: Grid convergence of the cascaded MRT LBM for the strain rate in a 2D Taylor-Green 
vortex flow with k = 1 and k = 2. 



19 



/ 



\ 
\ 
\ 
\ 
\ 
\ 
\- 
\ 
\ 
\ 
\ 
\ 
\ 



/ 



/ 



/ 



/ 



L 



/ 



X 



/ 



/ 



/ 



/ 



/////// A ////// / 



FIG. 4: Illustration of the geometry of a lid-driven cavity flow. 



test for CFD methods and has been a subject of many investigations using a variety of meth- 
ods (see e.g. [311435] ). Grid convergence for this problem has been studied using different 
collision models (SRT and standard MRT) for the LBM by various researchers (e.g. [12] )■ 
In this section, the aim is to analyze the grid convergence and an estimation of the order 
of accuracy of the cascaded MRT LBM for this flow problem. More detailed accuracy in- 
vestigation of the various flow features will be carried out in the next section. While the 
geometry is simple from the boundary condition implementation point of view, the flow 
contains singular points and becomes very complicated in terms of flow structures, partic- 
ularly as the Reynolds number increases (see e.g. [51] for a review). A schematic of the 
arrangement of the boundaries in a 2D lid-cavity flow is shown in Fig. |4} Fluid is enclosed 
inside a square cavity of length, L, and is set into motion by the moving upper wall that 
has a constant velocity Uo- The side and the bottom walls are considered to be stationary, 
which allows to implement a simple half-way bounce-back boundary condition on them. 
However, because the upper wall is in constant motion, a momentum correction needs to 
be added [SB] into the regular bounce-back scheme for the upper boundary. This is imple- 
mented as fa{i, Ny — 1) = fail, Ny — 1) + GpWaCayUp, whcrc fa{i, Ny — 1) is the post-colhsion 
distribution function, for a = 4,7,8, with a = 2,5,6 as the opposite directions of a, and 
Wa is the weighting factor [36]. Ghia et al. [31] have systematically studied this problem in 
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much detail by employing a vorticity-stream function formulation of the 2D incompressible 
Navier-Stokes equations, which is solved by a multigrid method. Some of their numerical 
results have been used for making accuracy comparisons in this work which will be discussed 
in a later section. Because of the lack of analytical solutions, the computed solutions ob- 
tained by a relatively very fine grid resolution, i.e. with i.e. 801 x 801, are treated as the 
approximate benchmark or reference ("analytical") solutions. Not only is the convergence 
of velocity fields tested, but also the grid convergence of the components of the strain rate 
tensor is considered. It may be noted that the study involving the latter quantity has not 
so far received enough attention for this problem using the LBM. 

The components of the velocity field and the strain rate tensor at the centerlines of the 
cavity in both vertical and horizontal directions are computed for a given Reynolds number 
once the solutions converge to steady state. The solutions are considered to reach steady 
state convergence when the relative global errors is small than 10~^^. Again, diffusive scaling 
is employed to set the parameters for different grid resolutions consisting of 13 x 13, 19 x 19, 
25 X 25, 31 X 31, 37 x 37, 49 x 49, 61 x 61, 85 x 85, 97 x 97 and 121 x 121 nodes. Figure [5] 
shows the grid convergence of the U-component of the velocity field at a Reynolds number 
of 100. It is found that the best fit slopes are 2.11 and 2.19 along the vertical and the 
horizontal centerlines, respectively, for the U-velocity. Likewise, the slopes are 2.18 and 2.11 
respectively along the vertical and the horizontal centerlines for the V-velocity as shown in 
Fig. 6 For the normal strain rate tensor component |^ , the slopes are found to be 1.81 and 
1.95 respectively for the vertical and the horizontal centerlines, which is shown in Fig. [7} 
Furthermore, it is seen that along the vertical and horiztonal centerlines, the strain rate 
tensor component |^ + ff has the slopes of 2.12 and 2.07, respectively, for grid convergence 
(see Fig. [s]). One reason why the slopes are either somewhat higher or lower than 2, rather 
than very close to the ideal value as seen with the other two problems discussed before, is 
that the reference solution for obtaining the relative error is taken to be that of the numerical 
solution with the very fine grid. This is often the practice as the "analytical" solution does 
not exist for this problem. 

Overall, it is seen that the cascaded MRT LBM gives a very respectable second order 
accuracy for a variety of fiows, including the relatively simple Poiseuille fiow and decaying 
Taylor-Green vortex fiow, and for relatively complex fiows such as the lid-driven cavity fiow. 
The method is found to be second order accurate not only for the velocity field, but also for 
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FIG. 5: Grid convergence of the cascaded MRT LBM for the U- velocity component in a 2D Ud- 
driven cavity flow for Re = 100. 




FIG. 6: Grid convergence of the cascaded MRT LBM for the V-velocity component in a 2D lid- 
driven cavity flow for Re = 100. 
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FIG. 7: Grid convergence of the cascaded MRT LBM for the strain rate tensor component ^ in 
a 2D Ud-driven cavity flow for Re = 100. 

their derivatives for the above problems. 

IV. ACCURACY STUDIES ON THE BENCHMARK PROBLEMS 

Let us now make a more detailed comparison of the accuracy of the solutions computed 
using the cascaded LBM with prior results involving either analytical or other numerical 
solution for the flow fields of the three benchmark problems considered in the previous 
section. 

A. 2D Poiseuille Flow 

Figure |9] shows a comparison of the velocity profiles of the 2D Poiseuille flow between 
the results obtained using the cascaded LBM and the parabolic analytical solution at a 
constant Reynolds number of 200 with constant relaxation time r = 0.515 for different grid 
resolutions in the wall normal direction starting from 26 to 401. Here, diffusive scaling is 
employed in the selection of parameters. That is, as the resolution is doubled, the maximum 
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FIG. 8: Grid convergence of the Ccisca^cied MRT LBM for the strEiin rate tensor component -1- -q^ 
in a 2D hd-driven cavity flow for Re = 100. 

flow velocity or the Mach number is decreased by a factor of 2. The results are in excellent 
agreement with the analytical solution, in which the maximum relative error is less than 
0.22 percent. 



B. 2D Decaying Taylor-Green Vortex Flow^ 



Using the same set of parameters specified for this time-dependent problem in the pre- 
vious section, we now compare the computed U— and V— velocity components along the 
vertical and horizontal centerlines, respectively, with the corresponding analytical solutions 



(Eq. (38)-(39)) at three different representative instants. Figures 10 and 11 show such a 
comparison of the velocity components at times t = 6.55T, 13.10T and 25.20T, where the 
characteristic time T is defined in the previous section, refiecting the decaying of the initial 
vortex distribution. It is evident that the cascaded MRT LBM is in excellent agreement 
with the analytical solution at all times shown. 
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FIG. 9: Comparison of the velocity profiles in a 2D Poiseuille flow for Re = 200 at different grid 
resolutions N in the wall normal direction with a constant relaxation time r = 0.55. 

C. 2D Lid-driven Cavity Flow 

Let us now consider more detailed features of the lid-driven cavity flow problem discussed 
in the last section at various Reynolds numbers in order to make quantitative comparisons. 



Figures 12 and 13 show the U- and V- components of the velocity, respectively, along the 
centerlines of the square cavity at Reynolds numbers of 100, 400, 1000, 3200, 5000, and 7500 
obtained using the cascaded MRT LBM along with the previous numerical data presented by 
Ghia et al [31]. The cascaded MRT LBM results corresponding to the finest grid considered 
earlier, i.e. for the 401 x 401 grid resolution are chosen to make comparison. In these figures, 
the solid lines represent the computed results obtained by the cascaded MRT LBM, and the 
symbols are the prior data provided by Ghia et al [31J. The velocities are normalized by the 
lid velocity Uq- Very good agreement is seen for all the Reynolds numbers considered. 

In a previous work, it was established that the standard MRT LBM based on raw mo- 
ments is superior when compared with the SRT LBM for the computation of lid-driven 
cavity flow [12j. Hence, it would be sufficient to make a direct comparison between the 
cascaded MRT LBM based on central moments and the standard MRT LBM for various 
flow characteristics of this problem. First, in order to provide a global characteristics of 
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FIG. 11: Comparison of the V- velocity component in a decaying Taylor-Green vortex flow for 
Re = 14.4 at three different non-dimensional times T: t = 6.55T, 13.10T and 26.20T. 
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FIG. 12: Comparison of the U- component of the velocity field along the vertical centerline of 
the cavity flow at various Reynolds numbers: Re = 100, 400, 1000, 3200, 5000 and 7500. Lines - 
cascaded MRT LBM and symbols - data by Ghia et al |31| . 




FIG. 13: Comparison of the V- component of the velocity field along the horizontal centerline of 
the cavity flow at various Reynolds numbers: Re = 100, 400, 1000, 3200, 5000 and 7500. Lines - 
cascaded MRT LBM and symbols - data by Ghia et al |31| . 
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the flow fleld, it would be interesting to compare the streamhnes in the cavity at various 
Reynolds numbers. It is known that at a certain Reynolds number above 7500, the flow 
field becomes unsteady and we restrict such comparisons for stationary state solutions only. 



Hence, Fig. 14 shows the computed streamlines at Reynolds numbers of 100, 400, 1000, 
5000 and 7500 using both the above methods. The streamlines computed by both these 
approaches are plotted side-by-side for comparison. It is found that the streamlines ap- 
pear to be remarkably very similar with both the raw moment and central moment based 
approaches. At Reynolds numbers of 100, 400 and 1000, a major vortex appears around 
the geometric center of the cavity with two minor vortices around the lower corners. Since 
the lid is driven from left to right, the major vortex circulates in a clockwise direction and 
the two minor vortices circulate in a counter-clockwise direction. At Reynolds numbers 
of 3200 and 5000, in addition to the vortices that exist with the lower Reynolds number 
cases, there appears another minor vortex on the left upper corner, which circulates in a 
counter-clockwise direction. When the Reynolds number increases further to 7500, a fourth 
minor vortex is found on the right lower corner, which circulates in a clockwise direction. 
All the above flow features correspond to steady states. Furthermore, in order to provide a 
more detailed comparison, we present various secondary vortices that appear in the cavity 



at Re = 7500 in Fig. 15 Again, remarkable similarity between the cascaded MRT LBM and 
the standard MRT LBM is found for these more detailed secondary flow structures. 

In order to provide a more quantitative perspective. Fig. [16] illustrates a comparison of 
the center of the primary vortex location in the cavity flow at different Reynolds numbers 
{Re = 100,400, 1000,3200,5000, and 7500) between the cascaded and standard MRT LBM 
as well as the data by Ghia et al [31]. From the earlier streamline plots, it can be observed 
that the location of the primary vortex moves towards the geometric center of the cavity 
as the Reynolds number increases. The computed results using the cascaded MRT LBM 
and the standard MRT LBM are in excellent agreement (within 0.014 percent) with each 
other for all Reynolds numbers. In addition, they are both in very good agreement with the 
data by Ghia et al [31j to within 0.50 percent for all Reynolds numbers. These quantitative 
results for the primary vortex locations are enumerated in Table |T] In addition. Table [TT] 
presents a comparison between the above two methods and the prior numerical data for 
the location of secondary vortices at different Reynolds numbers. Again, both the cascaded 
MRT LBM and the standard MRT LBM are in excellent quantitative agreement for the 
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(a) Cascaded MRT Re = 100 



(b) Standard MRT Re = 100 






(e)Cascaded MRT Re = 1000 



(f)Standard MRT Re = 1000 



location of these detailed secondary vortical structures with the data by Ghia et al |31] . 
Another useful global characteristic for comparison is the vorticity contours in the cavity 



at different Reynolds numbers. Figure 17 shows the vorticity contours computed using both 
the standard MRT LBM and the cascaded MRT LBM at three different Reynolds numbers 
{Re = 100,400, and 1000). As Reynolds number increases, the vorticity contours become 
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(g) Cascaded MRT Re = 3200 



(h) Standard MRT Re = 3200 




(i) Cascaded MRT Re = 5000 (j) Standard MRT Re = 5000 




(k) Cascaded MRT Re = 7500 (1) Standard MRT Re = 7500 

FIG. 14: Comparison of the streamlines in a 2D lid-driven cavity flow at diff'erent Reynolds numbers 
computed with cascaded (central moment) MRT LBM and standard (raw moment) MRT LBM: 
Re = 100, 400, 1000, 3200, 5000 and 7500. Solutions obtained using 201^ grids with both methods. 

denser and denser approaching the boundary walls. Overall, the vorticity distribution is 
found to be very similar using both the methods for all the Reynolds numbers considered 
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(a) Cascaded MRT Top 



(b) Standard MRT Top 





(c) Cascaded MRT Bottomleft 



(d)Standard MRT Bottomleft 





(e) Cascaded MRT Bottomright (f) Standard MRT Top 

FIG. 15: Comparison of the streamlines of the secondary vortices in a 2D lid-driven cavity flow 
at Re = 7500 computed with cascaded (central moment) MRT LBM and standard (raw moment) 
MRT LBM. 
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FIG. 16: Comparison of the Cartesian coordinates of the location of the center of the primary 
vortex in a Hd-driven cavity flow at different Reynolds numbers. 



TABLE I: Comparison of the location of the primary vortex in a lid-driven cavity flow at different 
Reynolds numbers. 



Re 


Cascaded MRT LBM 


Standard MRT LBM 


Ghia et a/(1982) [H 


100 


(0.61482,0.73543) 


(0.61467,0.73524) 


(0.61720,0.73440) 


400 


(0.55380,0.60514) 


(0.55380,0.60514) 


(0.55470,0.60550) 


1000 


(0.53070,0.56512) 


(0.53070,0.56512) 


(0.53130,0.56250) 


3200 


(0.51778,0.54027) 


(0.51777,0.54028) 


(0.51650,0.54690) 


5000 


(0.51499,0.53522) 


(0.51497,0.53524) 


(0.51150,0.53520) 


7500 


(0.51299,0.53186) 


(0.51298,0.53188) 


(0.51170,0.53220) 



thus corraborating the earher results. 

As discussed earher, one of the useful features of kinetic schemes such as the cascaded 
MRT LBM is that the components of the strain rate tensor can be obtained locally from the 



components of the non-equilibrium moments of the distribution function (see Eqs. (31 )-(33 )) 
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TABLE II: Comparison of the location of various secondary vortices in a lid-driven cavity flow at 
differnt Reynolds numbers. 



First Secondary Vortex 





Re 


Cascaded MRT LBM 


Standard MRT LBM 


Ghia et ai(1982) [ 




100 


NA 


NA 


NA 






NA 


NA 


NA 


Top 


1 non 

Qonn 
ozUU 


NA 

t^U.Uu4 ( , u.oy iv) 


NA 

l^u.uo4D, u.oy i o) 


NA 

l^U.U04; , U.oyo4J 




5000 


(n 0644 qosi 1 


(0 0641 9076') 


(0 0625 9102") 




7500 


(n 0676 9102") 


Co 0677 9099') 


Co 0664 9141") 




100 


(0.0387,0.0387) 


(0.0373,0.0373) 


(0.0313,0.0391) 






^^u.uuoo, u.u4:yo J 






Bottom Left 


1 non 

•^900 
OZUU 


(n ns49 n n7Qi ^ 

l^U. UO Z i , U. iZU / ) 


(f\ nsA9 n 0701 ^ 
/'o ns9i n 1 907^ 






5000 


(0.0740,0.1378) 


(0.0740,0.1378) 


(0.0703,0.1367) 




7500 


(0.0654,0.1536) 


(0.0654,0.1536) 


(0.0645,0.1504) 




100 


(0.9383,0.0658) 


(0.9386,0.0654) 


(0.9453,0.0625) 




400 


(0.8833,0.1243) 


(0.883,0.1243) 


(0.8906,0.1250) 


Bottom Right 


1000 
3200 


(0.8631,0.1128) 
(0.8229,0.0853) 


(0.8631,0.1128) 
(0.8229,0.0852) 


(0.8594,0.1094) 
(0.8125,0.0859) 




5000 


(0.8037,0.0739) 


(0.8037,0.0739) 


(0.8086,0.0742) 




7500 


(0.7892,0.0663) 


(0.7893,0.0663) 


(0.7813,0.0625) 


Second Secondary Vortex 




100 


NA 


NA 


NA 




400 


NA 


NA 


NA 


Bottom Left 


1000 
3200 


NA 

(0.0075,0.0075) 


NA 

(0.0073,0.0073) 


NA 

(0.0078,0.0078) 




5000 


(0.0075,0.0075) 


(0.0074,0.0074) 


(0.0117,0.0078) 




7500 


(0.0125,0.0125) 


(0.0115,0.0115) 


(0.0117,0.0117) 




100 


NA 


NA 


NA 




400 


(0.9926,0.0075) 


NA 


(0.9922,0.0078) 


Bottom Right 


1000 
3200 


(0.9923,0.0076) 
(0.9875,0.0113) 


(0.9928,0.0073) 
(0.9885,0.0115) 


(0.9922,0.0078) 
(0.9844,0.0078) 




5000 


(0.9775,0.0200) 


(0.9771,0.0193) 


(0.9805,0.0195) 




7500 


(0.9508,0.0429) 


(0.9509,0.0429) 


(0.9492,0.0430) 



Third Secondary Vortex 





100 


NA 


NA 


NA 




400 


NA 


NA 


NA 




1000 


NA 


NA 


NA 


Bottom Right 












3200 


NA 


NA 


NA 




5000 


NA 


NA 


NA 




7500 


(0.9964,0.0037) 


NA 


(0.9961,0.0039) 
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(a) Cascaded MRT Re = 100 



(b) Standard MRT Re = 100 




(c) Cascaded MRT Re = 400 (d) Standard MRT Re = 400 




(e)Cascaded MRT Re = 1000 (f)Standard MRT Re = 1000 

FIG. 17: Comparison of the vorticity contours in a 2D lid-driven cavity flow at different Reynolds 
numbers computed with cascaded (central moment) MRT LBM and standard (raw moment) MRT 
LBM: Re = 100, 400 and 1000. 
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The cavity flow being a shear driven problem generally has all the components of the strain 
rate tensor non-zero, and whose magnitudes can dramatically change with the Reynolds 
number. Hence, this problem provides a good test for the evalution of the accuracy of the 
computation of strain rate tensor by kinetic theory considerations, i.e. using non-equilibrium 



moments (Eqs. (|3l|)-(33 )). For the sake of comparison, we will make use of the standard 
second-order central differencing of the velocity fleld to obtain the usual direct estimation of 
the strain rate tensor components. In this regard, flow at two different Reynolds numbers are 
considered {Re = 100 and 1000) and the components of the strain rate tensor are obtained at 
flve different locations within the cavity using the above two methods, which are enumerated 



in Table |III[ As the Reynolds number is increased from 100 to 1000, the magnitudes of the 
strain rate tensor change signiflcantly, which are quite well captured by the kinetic approach. 
Indeed, remarkably the local computation using the non-equilibrium moments are in very 
good agreement with the flnite- difference estimation at various locations in the cavity for 
both the Reynolds numbers, with the maximum difference within 2 percent. This further 
demonstrates the numerical fldelity of the approach. In particular, such non-equilibrium 
moments based approach for the strain rate components can be used in the subgrid scale 
models for large eddy simulation of turbulent flows using the cascaded MRT LBM. 



V. NUMERICAL STABILITY STUDIES ON THE BENCHMARK PROBLEMS 

We will now discuss the results of numerical stability studies. Among the three benchmark 
problems discussed earlier, the lid-driven cavity flow presents the most stringent test since 
it is a fully 2D problem with boundaries containing singularity and the flow is shear driven. 
In fact, such a cavity flow problem was considered in detail to determine stability regimes 
of the SRT and the standard MRT collision models in a recent work [T^. Earlier, its three- 
dimensional counterpart was also considered from this viewpoint [37] . These studies have 
demonstrated the superiority of the use of multiple relaxation times in providing controlled 
additional numerical dissipation to enhance numerical stability on either coarser grids or at 
high Reynolds numbers when compared with the single relaxation time models. Hence, it 
is appropriate to consider the 2D lid-driven cavity flow to establish the stability regime of 
the cascaded MRT LBM in the context of other collision models. We now make a direct 
comparison of the maximum threshold Reynolds number for numerical stability of the SRT 
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TABLE III: Comparison of the components of the strain rate tensor computed using the local 



non-equilibrium moments (Eqs. (|31[)-(33)) and the finite-differencing (second-order central) of the 
velocity field with the cascaded MRT LBM at five different locations within the cavity for two 
different Reynolds numbers {Re = 100 and 1000). 



Re = 100 







Location 


Non-eqm. Moments 


Finite Difference 


Difference 








A 


(- 

\ 4 


§) 


2.416 X 10"* 


2.415 X 10-* 


0.044% 








B 


(- 
V 2 


f) 


1.711 X 10-5 


1.721 X 10-5 


0.613% 




dv 
dy 




C 


(- 
V 2 


§) 


1.850 X 10-"' 


1.848 X 10-* 


0.102% 






D 


^ 2 ' 


3L 

4 


-9.020 X IQ-^ 


-9.025 X 10-5 


0.057% 








E 


iSL 

4 


' 2 


-3.541 X 10"* 


-3.536 X 10-* 


0.125% 


Re = 100 








A 


(- 

\ 4 


2 '' 


3.516 X 10-5 


3.526 X 10-5 


0.300% 








B 


V 2 


4 


-4.344 X 10-* 


-4.342 X 10-* 


0.045% 


du 
dy 


+ 


dv 
dx 


C 


(- 
V 2 


2 


-3.368 X 10-* 


-3.363 X 10-* 


0.135% 






D 


(- 

2 ' 


3L N 

4 


4.599 X 10-* 


4.603 X 10-* 


0.093% 








E 


,3L 
4 


' 2 


-5.290 X 10-* 


-5.280 X 10-* 


0.198% 


Re = 1000 








A 


V 4 


f) 


4.220 X 10-5 


4.217 X 10-5 


0.050% 








B 


^ 2 


f) 


2.196 X 10-5 


2.209 X 10-5 


0.596% 




dv 
dy 




C 


V 2 


§) 


5.017 X 10-5 


5.017 X 10-5 


0.008% 






D 


2 ' 


31, \ 
4 > 


2.370 X 10-5 


2.372 X 10-5 


0.073% 








E 


/3L 

4 


-) 
' 2 


6.397 X 10-5 


6.446 X 10-5 


0.750% 


Re = 1000 








A 


V 4 


f) 


-1.344 X 10-* 


-1.334 X 10-* 


0.782% 








B 


\ 2 


f) 


8.137 X 10-5 


7.984 X 10-5 


1.914% 


du 
dy 


+ 


dv 
dx 


C 


V 2 


2 


-3.980 X 10-5 


-3.979 X 10-5 


0.021% 






D 


^ 2 ' 


3L\ 

4 


2.291 X 10-* 


2.289 X 10-* 


0.049% 








E 


4 


' 2 


-1.688 X 10-* 


-1.699 X 10-* 


0.626% 



LBM, the standard MRT LBM and the cascaded MRT LBM for this problem. With the 
cascaded MRT LBM, the relaxation parameters = UJ5 = 1/t are selected based on the 
specified kinematic viscosity, while the rest of relaxation parameters are set to unity for 
simplicity. For each approach, for a given grid resolution, the lid velocity was fixed and the 
relaxation time r was decreased gradually until the computation became unstable. 



Figure 18 shows the maximum Reynolds number {Re = UqL/v) that could be attained for 



each method before the computations became unstable, i.e. when the relative global error 
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increases rapidly or becomes exponentially large as the simulation progresses. Results are 
provided for different grid resolutions for these three approaches. It is clear that the cascaded 
MRT computations can reach Reynolds numbers that are about 2 or 3 times higher than that 
of the standard MRT approach and the standard MRT computations can reach Reynolds 
numbers that are 3 or 4 times higher than that of the SRT approach. The latter results are 
consistent with prior findings fi2\ [37] . Relaxation of different central moments at different 
rates provides a controlled additional numerical dissipation to maintain numerical stability. 
That is, maintaining frame invariance in conjunction with the use of multiple relaxation 
times further promotes the stability of the method. It may be noted that stabilization of 
certain classical methods have been achieved by constructing discretization operators that 
enforce Galilean invariance [T8H20]. Hence, it may be expected that explicitly incorporating 
an invariance property could aid with other standard mechanisms of stabilization of the 
LBM. As carried out in [12], we also perform an alternate stability test with the three 
approaches on a chosen coarse grid for this problem. In this test, the grid resolution is 
fixed at a relatively coarse resolution of 26 x 26, and then viscosity u (or equivalently r) 
is also set for all the three approaches. We then intend to find the maximum lid velocity 



which can maintain the stability of computations for 50,000 time steps [T2|. Figure 19 



shows how the three methods behave for this test. It is seen that the parameter regime 
or the maximum lid velocity for stability is considerably higher with the cascaded MRT 
LBM when compared with the other approaches. This further establishes the merits of 
the use of multiple relaxation times for central moment relaxation. Often, the stability of 
the CFD methods are characterized in terms of the grid or cell Reynolds number given 
by Rcc = UoAx/u (e.g. [38]). Thus, we also present the maximum cell Reynolds number 



for stability of the three approaches for this problem in Table |IV[ which demonstrates the 
advantages of the cascaded MRT LBM. 

Another important aspect is the computational cost. As shown previously, the cascaded 
MRT approach can be more stable with similar accuracy compared with the standard MRT 
for the lid-driven cavity flow. But if it is much more expensive for numerical computations 
than the standard MRT, its advantages will not be very useful. In this regard, we fully exploit 
all the optimization strategies that could be used with a moment approach, such as those 
specified in [SH] for the cascaded MRT LBM. It is found that for the 2D lid-driven cavity flow 
problem, the cascaded MRT LBM takes about 11.6% longer than the standard MRT LBM, 



37 







■ SRT 






■standard MRT 






^Cascaded MRT 


3 20000 


40000 


60000 80000 
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FIG. 18: Comparison of the maximum Reynolds number for numerical stability of different methods 
for simulation of the lid-driven cavity flow. 

TABLE IV: Comparison of the maximum cell Reynolds number {Rcc = UqAx/u) for numerical 
stability of different methods for simulation of the lid-driven cavity flow problem. 



Grid Resolution 


SRT LBM 


Standard MRT LBM 


Cascaded MRT LBM 


101 X 101 


14.14 


59.40 


148.50 


201 X 201 


14.21 


62.18 


165.83 


401 X 401 


14.25 


62.34 


199.50 



which is acceptable in view of the significant advantages in terms of numerical stability. It 
should be pointed out that these results pertain only to 2D problems. Additional work is 
required in three-dimensions to optimize the computational cost of the cascaded MRT LBM 
and also to optimize its relaxation parameters by means of a linear Fourier analysis. 
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1.8 1.84 1.88 1.92 1.96 2 

1/t 



FIG. 19: Alternative stability test to determine the maximum threshold lid velocity for different 
methods for a chosen coarse resolution (26 x 26). 

VI. SUMMARY AND CONCLUSIONS 

Galilean invariance is one of the main physical attributes in the description of the fluid 
motion. This is naturally achieved by considering dynamical changes in terms of central 
moments in kinetic schemes, as was done in the recently introduced cascaded LBM. Enforcing 
frame invariance is generally expected to have a positive influence on numerical stability as 
seen in some recent work with other classical schemes. The use of multiple relaxation times 
(MRT) in the central moment or cascaded LBM brings in the various flexibility associated 
with the standard MRT LBM based on raw moments. In particular, the relaxation of 
different central moments at different rates introduces additional dissipation as in the raw 
moment based approach, which can lead to enhanced stability. 

In this paper, we discussed our results from systematic numerical studies on grid conver- 
gence, accuracy, and stability of the cascaded MRT LBM. We have chosen three commonly 
used 2D benchmark problems including the Poiseuille flow, the decaying Taylor-Green vor- 
tex flow, and the lid-driven cavity flow. In the grid convergence tests, the cascaded MRT 
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approach has been found to be second order accurate under diffusive scahng for all the 
benchmark problems considered. These results are shown to hold not only for the velocity 
field, but also for the components of the strain rate tensors. Furthermore, comparisons of 
the numerical accuracy of the cascaded MRT LBM were made with other collision models 
and also with prior analytical or numerical results based on the solution of the Navier-Stokes 
equations. These demonstrated that the cascaded MRT LBM is in excellent agreement with 
the prior results for all the canonical problems considered. In particular, the detailed flow 
structures for the more complex lid-driven cavity flow predicted by the cascaded MRT LBM 
are in very good quantitative agreement with the standard MRT LBM. In addition, the 
utility and the accuracy of the use of non-equilibrium moments with the cascaded MRT 
LBM for the computation of the components of the strain rate tensor is demonstrated. Fi- 
nally, stability tests on a 2D lid-driven cavity flow problem was carried out, which showed 
substantial improvements in numerical stability of the cascaded MRT LBM, with higher 
threshold Reynolds numbers, when compared to other models. With the use of proper op- 
timization strategies, the 2D cascaded MRT LBM was found to be only about 10% to 20% 
more expensive when compared to the standard MRT LBM in terms of computational time. 

Future work could include further development of more optimized formulations of the 
three-dimensional cascaded LBM based on central moments with a view to maintain compu- 
tational efficiency and their applications to unsteady multiscale problems such as turbulence. 
Optimization of the relaxation parameters by a linear Fourier analysis to introduce adequate 
additional dissipation for enhanced numerical stability while maintaining necessary physics 
with this approach is also desired. 
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